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We provide a closed-form estimator based on the VARMA representation 
for the unrestricted multivariate GARCH(1,1). We show that all parameters 
can be derived using basic linear algebra tools. We show that the estimator is 
consistent and asymptotically normal distributed. Our results allow also to 
derive a closed form for the parameters in the context of tem poral aggregat ion 
of multivariate GARCH(1,1) by solving the equations as in Hafner 20081 ] . 
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1. Introduction 

Estimating a multivariate GARCH(1,1) model is a challenging task. The most common 
tool for this purpose is the quasi maximum likelihood (QML) estimator which requires 
rather sophisticated optimization techniques. In this paper we present a simple and fast 
method of moments which makes the estimation of the multivariate G ARCH (1,1) model 
more accessible. Our resul ts represent the multivariate generalization of the analytical 
results already achieved by Kristensen and Linton 2006] for the scalar case. 

Our estimator is consistent and, under additional assumptions on the moments, asymp- 
totically normal distributed. Due to the difficulties in estimating multivariate GARCH(1,1) 
models our estimator may then be used to provide a consistent initial estimate when im- 
plementing numerical optimization techniques for the QML estimation. This is especially 
true when large-scale models are employed. 

Several restricted models have been proposed by the previous literature in o rder to re- 
duce the number of parameters, suc h as Diagonal VEC ( Bollerslev et al. 1988]), BEKK- 



GARCH (jEngle and Kronerl |l995i |). CCC-GARCH dBollerslevI |l990^ . Interestingly 
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our results are valid in general. Therefore in the framework we stick to the unrestricted 
multivariate G ARCH (1,1). 



Finally, our results extend the results of lHafnerl 2008] in the context of temporal 
aggregation of multivariate GARCH(1,1). Indeed, our results allow to derive the pa- 
rameters of the temporally aggregated GARCH for any aggregation frequency. In other 
words, given the parameters of the disaggregated process, those of the aggregate one 
are analytical functions of the disaggregate parameters. Alternatively, one can also use 
the moments of the disaggregated GARCH to produce an initial estimate of the param- 
eters of temporally aggregated processes. The former estimator is again consistent and 
asymptotically normal when some moments conditions hold. 



moments 



2. Framework 

Consider the following unrestricted multivariate GARCH (1,1) model 

y t =Hl /2 e t , t = l,2,...,n, 

where yt is a d-dimensional zero-mean, serially uncorrelated process. In addition, we 
have that ej 6 M. dxl is an i.i.d. white noise vector with zero mean and variance Id- 
Moreover, the conditional covariance matrix is given by 



vech (H t ) =c + A\eck(yt-ivi-i) + B vech(^ t -i), t = 2, 3, 



,n, 



(1) garchrec 



where vech(M) represents the operator that stacks the elements of the lower triangular 
part of a symmetric matrix M to form a J X 1 vector, with d = Mi^hil, 
In what follows we make the following assumptions: 

1. Ht is positive definite almost surely for each t. 

2. All eigenvalues of the matrix A + B have modulus smaller than one. 

3. The process yt is ergodic, /3-mixing, and strictly stationary. 

4. The fourth moments of yt exist and are finite. 

Boussama 20061 ] provides sufficient conditions that ensure a strictly stationary, er- 



godic and /3-mixing solutio n of the vector GARCH process (these can also be found in 
Francq and Zakoi'an . 20ld . Theorem 11.5]). 



The following stronger assumption is used only in some central limit results: 
5 The eighth moments of yt exist and are finite. 



When the distribution of ej is spherical, Hafnerl . 120031 . Theorem 3] has given an algebraic 
condition equivalent to Assumption H] that is easy to test in practice. However, we do 
not need to assume sphericity here. 
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The VARMA(1,1) representation of the multivariate GARCH(1,1) is obtained by defin- 
ing ht = vech(H t ), xt = vech.(y t yj) and £t = %t — fh', the recurrence relation (JT]) is 
equivalent to 

x t = c + $x t -i+£t ~ BCt-i- 
By eliminating xt—i recursively, we find that asymptotically the following formula holds 

oo 

x t = h + ^2&i£ t -i, 

where <P = A + B; h = vech(H) = (I - $)- l c; <9 = I d and 9 t = (A + Bf^A for 
i > 1. The interest of this formulation lies in the fact that & is a martingale difference 
sequence. We define 17 = i£[£t£f]; note that E[£t£t+s] = for s > 1. 



sec : closedForm 



3. Closed form Estimation 

Under Assumptions 1-4, the autocovariances of x% exist and are finite, and they are 
given by 

oo 

M k = E[(x t+k - h)( Xt - h) T \ = e l+k sej. 

8=0 

From the VARMA(1,1) representation using the standard Yule- Walker results we have 
that 

M k+ i = <PM k , for all k > 1, 



thus tf> can be obtained analytically as 

& = M k+1 M k ~ 1 , 



for all k > 1. 



(2) 



Th ese results ar e well known and can also be found for example in the book of Reinsell 
1997] as well as in lHafnei [2003? ] page 32 (for the univariate case see lKristensen and Lintonl 
2006]). Therefore, we can estimate <P = M2M1 . Consider now the first-order moving 
average vector 

jt = x t - <Px t -i = c + it~ B^t-i- 



The autocovariances of jt are 

r = E[(j t - c)(j t - c) T ] = U + BUB T 



Mn 



Mi<Z> T - <f>Aff + <PM $ T , 



r 1= ,E[{jt-c){jt-i-c) T ] 



-BE = Mi - M 2 <P T - <PM + <2>Mi<2> T = Mi - <2>M . 



(3) 

We can combine the former two equations with simple manipulations to derive two 
separate equations for B and £ 



r{ + r B T + r 1 (B 
r 



T\2 



0, 



s + r x s- l rj. 



(4) 
(5) 



Phi 



Gammas 



pme 
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In t he scalar case, (d} is a quadratic equation; the approach suggested bv lKristensen and Linton 
2006] consists essentially in deriving an estimator by solving this equation. 



This method, however, need not be restricted to the univariate GARCH. In the multi- 
variate case, basic linear algebra techniques can be used to derive a closed form in terms 
of eigenvalues and eigenvectors. We present them in the next section. 

3.1. Closed formula for B 

The following procedure can be used to obtain B as analytical function of i~b an d A- 
1. Form the 2d x 2d matrix 



P 



/ 

_ r -i r r n-1 



(6) 



2. One can prove (Lemma in the following) that the eigenvalues of P come in pairs 
(A, 1/A). Therefore, unless there are eigenvalues that lie exactly on the unit circle, 
half of the 2d eigenvalues satisfy | A| < 1, and we may reorder them so that | Aj| < 1 
for i = 1, 2, . . . , d. Moreover, consider the associated eigenvectors Wi and partition 
them as 



IV; 



Ui.Vi e 



3. Now, a solution to the matrix equation (jl]) is given by 

"Ai 

B = {U T Y l DU T = L u 2 



A, 



A, 



Ui u 2 



(7) 



In Section and Appendix [Aj we recall the theoretical results in linear algebra that 
ensure the functioning of this procedure. 

3.2. Estimation using the closed formula 

Using this formula, an estimation procedure can be derived as following: 

1. Given the data, compute the observed average and the first three autocovariances 
of x t : 



1 



n r~f 



t=i 



1 n— 1 1 n— 2 

Mi = EK^+l - h)( x * ~ &) T ]> M 2 = t - h)(x t - hf). 

\ n ~ -v t=v in — 2) r~f 



1 n 

M^ = -Y)ixt-h){x t -h) T \, 



defP 
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Convergence 



2. Evaluate & = A^Mf 1 , f = M - Mi<L> T - Sm{ + &M S T , and A = M x 
as provided by ([3]). 



<PMn 



3. Use the above procedure based on eigenvalue computation to get an estimated B. 

4. Finally, recover the other two parameters as A = — B and c = (I — A — B)h. 

3.3. Asymptotic properties 

In this section and the following, the symbols A and ^ stand for convergence in prob- 
ability and distribution, respectively. 

The consistency of the QM L estimator has been shown by Jeantheau 19981 ] while 
Comte and Lieberman 20031 ] provides the asymptotic normalit y of the QML e stima - 



tor in the context of BEKK formulation. However, as noted in Bauwens et al. 20061 ]: 



The asymptotic properties of ML and QML estimators in multivariate GARCH mod- 
els are not yet firmly established, and are difficult to derive from low level assump- 
tions [...] Asymptotic normality of the QMLE is not established generally. [...] Re- 
searchers who use MGARCH models have generally proceeded as if asymptotic normality 
holds in all cases. Here we provide the asymptotic properties of our closed-form es- 
timator A = (c, vec(A), vec(-B)), which is function of the moments of x t only, that is, 
A = F(h,M ,M 1 ,M 2 ). 

Lemma 1. Let 

h 

vec Mq 
vec Mi 
vec M2 

Under Assumptions 1-4; rh A m. In addition, if Assumption 5 holds, we have 

dist 



m 



y/n(m 



m 



N(0,\P), 



(8) 



where & = Cov[m]. 



Proof. Under these assumptions, both h and Mk are finite, in addition h and Mk are 
consistent given the law of large numbers for stationary and ergodic processes. Finally, 
the joint distribution converges to a normal distribution thanks to the standard central 
limit theory for strongly mixing sequences. □ 

Note that m is consistent even if Assumption 5 does not hold; in this case, howeve r, 
m converges with a slower rate (see the discussion in Kristensen and Linton 2006]). 
Moreover, the limit joint distribution of the moments is not Gaussian . 

If the noise distribution is spherical, then, as noted by iHafnerl \2QQ$\ . the cross- 
covariance between h and vec(M/%) vanishes for each k, as it is an odd-power function of 
the noise et- In this case, the explicit expression for 'F is 



Cov[fc] 0, 
%d^d Cov[M fc ,M ; ] 



dx3d 2 
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Normality 



where: 
and 

Cav[M k ,Mj\ 



Cov[/i] =M + J2 M k + Jl M l 



k=l 



k=l 



E 



(vec [(x t+k - h)(x t - h) T ] - vec M h \ (vec [(x t+i - h)(x t - h) 



vec 



MA 



for k,l = 0,1,2. 

The following theorem represents a direct consequence of Lemma [TJ 

Theorem 2. Under Assumptions 1~4, the estimator A = (c,vec(A),vec(B)) = F(m) is 
consistent, i.e. A A A. In addition, if Assumption 5 holds, \/n(A — A)— >N(0, S) with 



( dF(m) 
\ dm 



/ dF(m) 
\ dm 



(9) 



Proof. The theorem follows from the continuous mapping theorem; in addition, the 
explicit expression for the covariance is obtained using the delta method. It remains to 
show that the partial derivative 9I q^ exists, which is proven in the end of Section □ 

Note that the asymptotic properties of the closed fo rm estimator might b e employed 
to prove these of the QML in general (as discussed by Bauwens et al. 20061 ] ). However, 
we leave this for future research. 



usesPartialDer 



4. Temporal aggregation 

An interesting consequence of the results above is the direct extension to the temporal 
aggregation. In fact, we can now derive (as well as estimate) the paramete rs of the 
temporally aggregated multivariate G ARCH (1,1) as discussed in Hafner 20081 ]. 



Temporal aggregation of a GARCH can be conducted in two different forms, depending 
on whether we are interested in stock or flow variables. We are interested in deriving a 
GARCH representation for the process y^ m ' aggregated over m periods, which is defined 
in the two cases as 

Vmt = Vmt, (stock variables) 

Vmt = Vmt + Vmt-i ■■■ + Vmt-m+1 + . (flow variables) 

We denote x^ = vech(y^ yl^ T ) . While in the stock case x^ = x m t, the relation 
between the second moments is more involved for flow variables. 



Hafnerl 2008] shows that the temporally aggregated process follows a weak VARMA(1,1) 
Tint +^ X m(t-1)+W D Sn(t-l)- u " 
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Then, equations for £>( m ) formally analogous to © are derived: 

r[ m) = -b^ s^ m \ 

with 



r (m) 
1 


m 




i=0 


7-. (m) 


— t s y 










7 s 


= -$ m ~ l B, 



for the case of stock variables, and 



r (m) 
1 


2m- 1 

8=0 
m— 1 

i=0 




r (m) 
1 1 






= k. 






= Ij+ A + $A -\ h # _1 A i = 


1, • • • , m — 1 


7 J 


= [ij + + ■ ■ ■ + (p m ~ 2 ]A - $ m - l B, 


i = m 




= + $i-m+l _| |_ $rn-2^ A _ 




J 2m-1 







for the case of flow variables (where an explicit expression for U™ can be found in 
Hafneri 2008 ] eq.19). The author, however, rearranges them to eliminate in a form 
that differs slightly from our @, and for which deriving an explicit solution is more 
complicated. In his words, "As of B^ m \ (29) is a system of nonlinear equations that 
cannot be solved explicitly.". With the too l s pro vided in this paper, an explicit solution 
is now available. Equation (29) in Hafner 20081 ] can be replaced with Q in this paper. 
Therefore using 



p(m) 



r (m)-l r (m)T 
1 1 1 1 



-r- 



i 

(m)— 1 j-,{ra) 
1 i 







we can carry on the procedure described in Section Again, the eigenvalues to choose 
are those inside the unit circle. 



defPm 
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alityAggregate 



rlafnerl [20081 ] shows that 2~o and r\" l> are analytical functions of A, B and E(£tC, 



-,(m) 



and these are are function of the moments of xt only. It turns out that a closed form 
estimator of the temporally aggregated G ARCH (1,1) can be derived as an analytical 
function A^ = (c^ m \ vec A^ m \vec B^) = G^ (h, vec M , vec M u vec M 2 ). In par- 
ticular, we can use for the estimation of the aggregated GARCH(1,1) the estimated 
moments of the high-frequency data, for which more information is available. 
The estimator enjoys the same asymptotic properties. 

Theorem 3. Under Assumptions 1~4, the estimator A^ = (c^ m \vec A^ m \ vec-E>( m )) = 
G( m )(m) is consistent, i.e. A^ A A^ m \ In addition, if Assumption 5 holds as well, we 
have: ^(A^ - A)^-N(0, sH), with 



'dG^{m)\ fdG^{m) 



dm 



dm 



(12) 



usesPartialDer 



Hafner and Romboutsl 2007 ] discuss the estimation of temporally aggregated multi- 



variate GARCH(1,1). Our estimator can be employed as a simple estimator when the 
number of observations is sufficiently high. Alternatively, it can be used as a consistent 
starting value for the QML estimation. 



: linear Algebra 



5. Palindromic matrix equations and eigenvalue problems 

In this section, we present the linear algebra results that lead to the estimator of B 
(as well as B^). These m atrix equations have b e en studied by many au t hors in linear 
algebra literature, see e.g. iGohberg et all |l982l |. lEngwerda etaH |l993l ]. iMeinil [20021 ] 
and the references therein; we focus here on providing a non-technical exposition. Proofs 
of some of the results are presented in Appendix lAl 

The problem of computing one or more pairs (A, u) satisfying 



solGenEig 



(a 2 a + xr + r{)u =o, 



(13) 



is known as palindromic quadratic eigenvalue problem iMackev et al.l [20061 ] . The complex 
numbers A are called generalized eigenvalues and the vectors u generalized eigenvectors. 
It is indeed a generalization of the standard eigenvalue problem, i.e., given a matrix A 
finding pairs (A, u) satisfying (A — \I)u = 0. 

First, we show that all the solutions to ([3D can be constructed from generalized eigen- 
values and eigenvectors of (fT5|) . 

Lemma 4. Let (Ai, Mi), (A2, u%), . . . , (Aj, u£) be d different pairs of generalized eigenval- 
ues and eigenvectors of the problem (|13p . such that the matrix U (as in ^) is invertible. 
Then, B = (C/ T )- 1 D[/ T is a solution of ©. All the solutions of © can be obtained in 
this way. 



pep 



Moreover, in our problem the possible values of A can be "paired". 
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pairinglemma Lemma 5. Let i~b o,nd A be real matrices, with i~b symmetric. If X ^ is a generalized 
eigenvalue of (|13|) . then 1/A is one as well. 



analytical 



Finally, the following result shows that we can reduce the palindromic eigenvalue 
problem to a standard eigenvalue problem. 



eigenlemma Lemma 6. Suppose that the matrix T\ is invertible. 



1. Let (A, u) be a solution of (|13|) , Then, A is an eigenvalue of P, with corresponding 
eigenvector 

u 



w 



Xu 



2. Conversely, if X is an eigenvalue of P with eigenvector w , then (A, u) satisfy (fT3"j) . 
where u is the vector formed by the first d entries of w. 

By combining all the above results, we can obtain a closed form for B. Note that 
different solutions are possible; namely, every choice of d eigenvalues out of the 2d of 
P gives a different B satisfying 0; however, only the one with |Aj| < 1, % = 1, 2, . . . , d 
results in a B with all its roots inside the unit circle. 

Remark 7. The invertibility of i~i is not a crucial assumption. If I\ is singular, we can 
obtain in a similar way not an eigenvalue problem of the form Mv = Xv, but the slightly 
more general form Mv = XNv. This is known as generalized eigenvalue problem, and 
there are plenty of algorithms to find a closed- form solution to it. For instance, the 
Matlab command eig(M,N). Similarly, if P is not diagonalizable, solutions B to the 
matrix equation @ can be defined in terms of its Jordan canonical form. 

The existence of the partial derivatives ^[ fe ' vec , A /° ' vec » J ^ 1 ' wec ^ 2 \ that are needed in The- 

^ o(ft,vec Aiojvcc Aii ,vcc M2) 

orem[2]can be shown, again under the condition that P has no unimodular eigenvalues. 

Lemma 8. Suppose that the matrix P has no eigenvalues on the unit circle. Then, B 
is an analytical function of the equation coefficients Tq, 

Since i~b and T± are in turn analytical functions of the Mj, the partial derivative 
exists. A sketch of proof of this result is in the appendix, together with a more explicit 
expression for the Jacobian. 



6. Small sample issues 

unimodular 

The results provided in this paper should be employed with caution whenever the sample 
size n is not large enough. The closed-form estimator is based on the sample estimates 
of Mq, Mi, Mi. However, it may be the case that the sample moments do not respect 
all the stated assumptions. More specifically, three different kind of issues can arise: 

Positivity c, A and B do not guarantee that Ht is positive definite. 

Stationarity The roots of 3 = Mk+\M^ 1 lie on or outside the unit circle. 
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Invertibility B has unimodular eigenvalues (i.e., on the unit circle). 



When the GARCH parameters are estimated via maximum likelihood, the constraints of 
respecting these conditions are usually imp osed when solving the optimization problem; 
see for instance IChretien and Ortegal 20121 ] . Since black-box optimization routines are 
used, additional constraint are easy to impose, but they make the resulting problem 
more complicated to solve. On the plus side, they guarantee that the resulting model 
has the desired properties, provided that the iterative optimization procedure does not 
fail. Instead, with an exact moment-based estimator, if one or more of these conditions 
fail, then the best way out is modifying the sample moments or the estimates a posteriori 
to make sure that they satisfy these constraints. We discuss briefly these problems that 
may arise when our closed form estimation is employed. 



6.1. Positivity 



Sufficient conditions for pos it ivity are discussed by severa l authors (see lGourierouxl [19971 ] . 
Chretien and Ortega 20121 ] . Francq and Zako'ian 201d( ]). However, as far as we know, 
the problem of finding necessary conditions has not been dealt with in literature. Indeed, 
even the simpler problem of finding all linear maps among symmetric matr ix spaces that 
preserve positive semi-definiteness has no simple solution, see for instance [Bhatial . 120071 . 
Chapters 2 and 3]. Our estimation procedure does not always produce an estimated 
GARCH satisfying the sufficient conditions cited above. We do not deal here with the 
problem of finding a weaker set of conditions that can be preserved. 



6.2. Stationarity 



In small samples the estimate of ^ = Mk+\M^~ can have eigenvalues on or outside 
the unit circle; moreover, the values computed by choosing different values of k in 
the former expression will in general be different. The choice described above of taking 
5>W and ignorin g all the other autocoyarianc es ratios is the simplest way out of the 
latter problem. iKristensen and Lintonl 2006] discuss this problem in the scalar case, 
and suggest as another valid approach taking | (#W -|- ^>( 2 ) -j-^C 3 )^ or j n general any 

convex combination J2 w i@^- 

To avoid problems with noninvertibility or outliers, in the multivariate case it is more 
advisable use instead a least-square solution of the system 



uiiM\ W2M2 



w\Mi W2M3 



w n -iM n 



again with suitably-chosen weights Wi. None of these solutions (and no choice of weights) 
clearly stands out. In particular, all of them may result in estimates with eigenvalues 
equal or larger than 1. When this happens, a simple fix is projecting the estimate on the 
space of acceptable GARCH solutions by altering the eigenvalues that lie on or outside 
the unit circle. 
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6.3. Invertibility 

Our assumptions on the solution B guarantee that it has all eigenvalues inside the unit 
disc, and thus that the matrix P has no unimodular eigenvalues. However, once again the 
sample autocovariances from a finite-time realization of a GARCH process may lead to 
sample estimates of Jo and P\ that do no necessarily guarantee that P has no eigenvalues 
on the unit circle. This is especially true when the process is close to a non-invertible 
one. When B has unimodular eigenvalues the invertibility condition does not hold. In 
addition the following theorem sheds light on the consequences of having eigenvalues 
lying on the unit disk. 



Theorem 9. The following results hold. 

1. If P has no eigenvalues lying on the unit circle, then there exist unique solutions 
B and E, where U = U T and |A| < 1 for each eigenvalue A of B, and they can be 
computed with the above procedure. 

2. If P has eigenvalues lying on the unit circle, then B and a positive definite £ satis- 
fying ([2]) exist only if some strong additional conditions are satisfied ( in particular, 
all unimodular eigenvalues should have even multiplicity). In this case, B always 
has unimodular eigenvalues. 

A full proof is more technical th an those for the other linear algebra results that we 



reported; we omit it and refer to lEngwerda et al.l 1993] for a complete presentation. 



However, the last assertion is clear in view of our derivation: since the eigenvalues of 
B are a subset of those of P, B cannot have all its eigenvalues inside the unit circle 
if P has less than d eigenvalues in that domain. If the existence conditions are not 
satisfied, we can still compute solutions B with p{B) = 1 with the procedure of Section^ 
there are multiple solutions, according to which unimodular eigenvalues we choose, but 
none of them will result in a symmetric 27. Ad-hoc modifications of P can be made 
when unimodular eigenvalues are detected, but in general the accuracy of the computed 
solution is expected to decrease. Indeed, we show in Appendix [B] that the derivative of 
B with respect to the moments can become unbounded when it has eigenvalues equal 
to 1. 

We are currently working on developing a general procedure for computing a small- 
norm modification of the that makes the estimated model invertible, rooted on results 
in linear algebra and eigenvalue perturbation theory. 



We point out that the same problem arises in the scalar case treated bv lKristensen and Linton 



2006]: when (in our notation) Jb/A > 2, the scalar quadratic equation ([3]) has two com- 



plex conjugate solutions with modulus 1, and the procedure breaks down. 
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secproof s 



A. Proofs 
Proof of Lemma [4] 

Proof. Using the fact that the (Aj, Uj) are generalized eigenvalues, one can check directly 
that each column of the matrix 

r x UD 2 + r UD + r?u 

is zero; therefore, 

o = (ac/d 2 + r UD + rfu)u- x = riUDU-^iUDU^ 1 ) + r^uDU- 1 ) + if, 

as required. For the converse implication, let B T = (UDU^ 1 ) be the spectral decom- 
position of a solution; we can reverse all the steps and obtain that each (Aj,-Uj) is a 
generalized eigenpair. □ 

Proof of Lemma [5] 

Proof. Let A satisfy (Q2J) for some choice of u ^ 0. Since i~b, A are real, we can take the 
complex conjugate of every term and get 

(A 2 A + AT + rf)u = 0, 

where A and u denote (componentwise) complex conjugation. In particular, this implies 
that 

det(A 2 A + A\Tq + if) = 0. 
Then the determinant of its conjugate transpose must be as well, and thus 

det(A 2 A T + AA> + A) = 0. 



Multiply everything by p-, to obtain 



= det ^ (A 2 A T + \r + A) = det (if + jl + * I^j . 

Since this determinant is zero, the matrix is singular and there must be a vector u such 
that 

/ i /i\2 1 

-iT 1 - I 1 



But this equation shows that the pair fx>^) * s a ^ so a generalized eigenpair of the poly- 
nomial eigenvalue problem. □ 
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Proof of Lemma [5] 

Proof. 1. Let A,ti be a solution to (|13|) . We may verify explicitly that the equation 

/ 

holds. 

2. Let w be an eigenvector of P with eigenvalue A, and partition it as 





u 


= A 


n 




An 


An 



w 



From the eigenvalue condition Pw = Xw we obtain 

/ 

Expanding the two blocks and eliminating v from the resulting equations, one gets 





u 


= A 


u 










V 




V 



(USD- 
Proof of Lemma M 



□ 



The proof follo ws from some c l assica l results in matrix polynomials that can be found, 
for instance, in iGohberg et al.l 19821 ] . We give the sketch of a self-contained proof here. 
We start from a classical result in complex analysis, the Cauchy integral formula 



1 

27d 



1*1=1 



(z - A)" 1 dz 



'1 |A|<1, 
|A|>1. 



From this, a matrix version of the same integral follows for diagonal matrices 

1 



(zl - D)- 1 dz = 77, 

where II is the diagonal matrix such that Ha is zero if \Da\ > 1 and one if \Da\ < 1. Now 
a change of bases in both sides of the equation yields for all diagonalizable A without 
unimodular eigenvalues 

Ha, (14) 



— / {zI-A^dz 
2m J\z\=i 

with Ua the projector on the invariant subspace of A associated to the eigenvalues inside 
the unit circle. 

We may generalize further this formula to all A without unimodular eigenvalues, re- 
moving the diagonalizability of A from the requirements. Indeed, for a non-diagonalizable 
A, let us consider a sequence of matrices A^, each of them diagonalizable, that converge 
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uniformly to A. Such a sequence exists because diagonalizable matrices are dense in 
C nxn . Since the integrand function is olomorphic on the integration contour, limits 
and derivatives can be moved inside the integral. This shows that II a is an analytical 
function of A. 

In particular, we apply the formula for A equal to the 2d x 2d matrix P in (0), for 
which lip has rank d due to the eigenvalue pairing. Let us take any 2d x d matrix V 

such that IIpV has full rank and spans the range im lip. Since 



for the same subspace, it follows that 



B 1 



is another basis 



h o 



UpV 



(15) 



BTf ormula 



Since invertibility is a condition that holds on an open domain, (|15p holds locally with 
a constant V and provides an analytical expression for B in terms of Bp and thus of P. 
The matrix P is in turn a function of Zq> A- 



appDer 



B. Expressions for the derivatives 

In this section, we give a computable form for the Jacobian of F : (h, Mq, Mi, M2) 
(c, A, B), the function considered in Theorem [21 Rather than using vectorization to give 
an unwieldy matrix expression, we focus on describing its action as a linear map that 
takes a first-order perturbation of the moments (denoted by (h, Mq, M\, M2)) to one of 
the parameters (c,A,B). We shall use several times the expression for the derivative of 
the matrix inverse ^(M -1 ) = -M -1 (^M)M -1 . 

The relation between (h, Mo, Mi, M2) and A, A is easy to compute, by simply differ- 
entiating 0: 



A 
A 



Mi 
M 



<PM - <2>M , 



Mi <2> J 



Mi<F - *Aff - <2>Mf + <t>M $ T + <2M <Z> T + <2>M <F, 



(16) 



derl 



with ^> = M 2 Mf 1 - M 2 M^ 1 M 1 M^ 1 = M 2 Mf 1 - ^MiMf 1 , obtained by differentiat- 
ing d2D for A; = 1. 

We now differentiate (0) to obtain 



t - beb t = t- f x e- x es- x f1 = f - rxS^rf - riE- 1 ^ = i + /i# T + Btf. 

^17) 

This is a discrete-time Lyapunov equation (see for instance iGaiic and Qureshil 199a ] ) 
for U, which can be solved in closed form by vectorization or numerically by procedures 
such as Matlab's dlyap. The equation is uniquely solvable since we are assuming that 
P (B)<1. / ' _ 

Once we have E, we differentiate B = —TiS to obtain 



der2 



B 



-A^ 1 + FlS^SIT 1 = -(A + BE)!]' 1 . 



(18) 
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The derivatives of the remaining two parameters are given by A = <P — B and c = 

(J - 4)h + (i - &)L 

Putting together (fTUj) . CL7D, (USD, one can get to an expression for (c, A, B) as a function 
of (h, Mo, M%, M2). It does not look like there are any significant simplifications in the 
resulting expressions. The main message to infer from this computation is that the norm 
of (I — B® B)^ 1 , which appears when solving the discrete-time Lyapunov equation, has 
an impact on the magnitude of the derivatives; the closer B is to having unimodular 
eigenvalues, the more ill-conditioned the solution becomes. 
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